##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(ggpubr)
library(patchwork)
library("scales")

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--sample_list_public_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--maf_public_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"

    sample_list_file <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    sample_list_public_file <- paste(work_dir,"/public_ref/combine/MutationInfo.combine.tsv",sep="")

    rsem_file <- "~/20220915_gastric_multiple/dna_combinePublic/mRNA/CombineTMM.DNAUse.NJMU_TCGA.tsv"
    maf_public_file <- paste(work_dir,"/maf_public/All_use.addVAF.maf",sep="")
    
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/expression"
    gene <- "MUC6"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
out_path <- opt$out_path
rsem_file <- opt$rsem_file
sample_list_public_file <- opt$sample_list_public_file
maf_public_file <- opt$maf_public_file
gene <- opt$gene

dir.create( out_path , recursive = T )

##########################################################################################

info <- data.frame(fread(sample_list_file))
info_public <- data.frame(fread(sample_list_public_file))
dat_expression <- data.frame(fread(rsem_file))
dat_maf_public <- data.frame(fread( maf_public_file ))

##########################################################################################

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

check_gene <- c("DSC2" , "DSG2")
gene <- "MUC6"

dat_expression <- subset( dat_expression , gene_id %in% check_gene)
dat_maf_public <- subset( dat_maf_public , Hugo_Symbol == gene & Variant_Classification %in% Variant_Types )
colnames(dat_expression) <- gsub( "[.]" , "-" , colnames(dat_expression) )

##########################################################################################

info_public <- subset( info_public , From != "NJMU" )
info_public$ID <- info_public$Tumor
info_public$Class_sub <- info_public$Class
info_use <- rbind( info_public[,c( "ID" , "Tumor" , "Class" , "Class_sub" , "From")] , info[,c( "ID" , "Tumor" , "Class" , "Class_sub" , "From")] )

mutTumor <- unique(dat_maf_public$Tumor_Sample_Barcode)
info_mut <- subset( info_use , Tumor %in% mutTumor )
info_mut <- paste0(info_mut$ID , "_" , info_mut$Class_sub)
info_mut <- info_mut[info_mut %in% colnames(dat_expression)]

info_wild <- subset( info_use , !(Tumor %in% mutTumor) )
info_wild <- paste0(info_wild$ID , "_" , info_wild$Class_sub)
info_wild <- info_wild[info_wild %in% colnames(dat_expression)]

##########################################################################################
plotTpm <- function(tmm_combine = tmm_combine , out_name = out_name , title = title , width = width , y_tmm = y_tmm){

    my_comparisons_1 <- list( c(1, 2) )

    y_max <- y_tmm
    if(y_max > 2000){

        if(y_max > 10000){
            y_breaks <- 10000
        }else if(y_max > 5000){
            y_breaks <- 5000
        }else{
            y_breaks <- 2000
        }
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.1) , breaks = c( 0 , 500 , 1000 , 2000 , y_breaks ) , trans = sqrt_trans() )
    }else{
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.1))
    }

    plot <- ggplot( tmm_combine , aes( x = Type , y = TMM , color = Type ) ) +
        geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
        scale_fill_npg()+
        scale_color_npg()+
        facet_grid(.~Class)+
        xlab(NULL) +
        ylab(paste0("TMM of " , gene))+
        theme_bw() +
        ggtitle(title) +
        y_lab_lim +
        stat_compare_means(comparisons = my_comparisons_1) +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold',hjust=0.5,vjust=0.5),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.text.x = element_text(size = 10,color="black",face='bold',angle = 45 , hjust=0.5,vjust=0.5) ,
            axis.line = element_line(size = 0.5)) 
    return(plot)
}

getTMM <- function(exp_use = exp_use , type = type){
    sample <- sapply( strsplit(names(exp_use) , "_") , "[" , 1)
    class <- sapply( strsplit(names(exp_use) , "_") , "[" , 2)
    class <- sapply( strsplit(class , "-") , "[" , 1)
    tmm <- as.numeric(exp_use)
    tmp_dat <- data.frame( Sample = sample , Class = class , TMM = tmm )
    ## 合并一个人的同一样本
    tmp_dat <- tmp_dat %>%
    group_by( Sample , Class ) %>%
    summarize( TMM = median(TMM) )
    tmp_dat <- data.frame(tmp_dat)
    tmp_dat$Type <- type

    return(tmp_dat)
}

describeGene <- function(dat_expression = dat_expression , image_path = image_path ){

    tmp_exp <- dat_expression
    ## 突变型样本的表达
    type <- "MUC6_Mut"
    use_sample <- info_mut
    exp_use <- tmp_exp[use_sample]
    tmm_mut <- getTMM(exp_use = exp_use , type = type)

    ## 野生型样本的表达
    type <- "MUC6_Wild"
    use_sample <- info_wild
    exp_use <- tmp_exp[use_sample]
    tmm_wild <- getTMM(exp_use = exp_use , type = type)

    ## 合并
    tmm_combine <- rbind( tmm_wild , tmm_mut )
    tmm_combine$Class <- factor( tmm_combine$Class , levels = c( "IM" , "IGC" , "DGC" ) , order = T )

    y_tmm <- max(tmm_combine$TMM) * 1.1

    from <- "All"
    tmp_dat_use <- tmm_combine
    width <- 5
    title <- paste0(
            "Gene : " , gene , "\n" , from
    )
    out_name <- paste0( image_path , "/" , gene , ".MutVsWild.",from,".pdf" )
    plot <- plotTpm(tmm_combine = tmp_dat_use , out_name = out_name , title = title , width = width , y_tmm = y_tmm)
    ggsave(file=out_name,plot=plot,width=width,height=5)
    out_name <- paste0( image_path , "/" , gene , ".MutVsWild.",from,".tsv" )
    write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)

    ## TCGA和NJMU的图放一起
    width <- 5
    from <- "NJMU"
    index <- grep( "TCGA" , tmm_combine$Sample , invert = T )
    tmp_dat_use <- subset( tmm_combine[index,] , Class!="IM")
    title <- paste0(
            "Gene : " , gene , "\n" , from
    )
    plot1 <- plotTpm(tmm_combine = tmp_dat_use , out_name = out_name , title = title , width = width , y_tmm = y_tmm)
    out_name <- paste0( image_path , "/" , gene , ".MutVsWild.",from,".tsv" )
    write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)

    from <- "TCGA"
    index <- grep( "TCGA" , tmm_combine$Sample)
    tmp_dat_use <- subset( tmm_combine[index,] , Class!="IM")
    title <- paste0(
            "Gene : " , gene , "\n" , from
    )
    plot2 <- plotTpm(tmm_combine = tmp_dat_use , out_name = out_name , title = title , width = width , y_tmm = y_tmm)
    out_name <- paste0( image_path , "/" , gene , ".MutVsWild.",from,".tsv" )
    write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)

    out_name <- paste0( image_path , "/" , gene , ".MutVsWild.From.pdf" )
    ggsave(file=out_name,plot=plot1 + plot2,width=width,height=4)
}

image_path <- out_path

for( gene in unique(dat_expression$gene_id) ) {
    dat_expression_use <- subset( dat_expression , gene_id == gene )
    print(gene)
    describeGene(dat_expression = dat_expression_use , image_path = image_path )
}

##########################################################################################
## 只看肠化里面
tmm_combine_final <- c()
for( gene in unique(dat_expression$gene_id) ) {
    tmp_exp <- subset( dat_expression , gene_id == gene )
    
    ## 突变型样本的表达
    type <- "MUC6_Mut"
    use_sample <- info_mut
    exp_use <- tmp_exp[use_sample]
    tmm_mut <- getTMM(exp_use = exp_use , type = type)

    ## 野生型样本的表达
    type <- "MUC6_Wild"
    use_sample <- info_wild
    exp_use <- tmp_exp[use_sample]
    tmm_wild <- getTMM(exp_use = exp_use , type = type)

    ## 合并
    tmm_combine <- rbind( tmm_wild , tmm_mut )
    tmm_combine$Class <- factor( tmm_combine$Class , levels = c( "IM" , "IGC" , "DGC" ) , order = T )
    tmm_combine$gene <- gene
    tmm_combine_final <- rbind( tmm_combine_final , tmm_combine )
}

tmm_combine_final <- subset( tmm_combine_final , Class == "IM" )
tmm_combine <- tmm_combine_final

y_tmm <- max(tmm_combine$TMM) * 1.1
tmp_dat_use <- tmm_combine
width <- 5

my_comparisons_1 <- list( c(1, 2) )

y_max <- y_tmm
if(y_max > 2000){

    if(y_max > 10000){
        y_breaks <- 10000
    }else if(y_max > 5000){
        y_breaks <- 5000
    }else{
        y_breaks <- 2000
    }
    y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.1) , breaks = c( 0 , 500 , 1000 , 2000 , y_breaks ) , trans = sqrt_trans() )
}else{
    y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_tmm*1.1))
}

plot <- ggplot( tmm_combine , aes( x = Type , y = TMM , color = Type ) ) +
    geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6) +
    geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
    scale_fill_npg()+
    scale_color_npg()+
    facet_grid(.~gene)+
    xlab(NULL) +
    ylab("TMM")+
    theme_bw() +
    y_lab_lim +
    stat_compare_means(comparisons = my_comparisons_1) +
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='none',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        plot.title = element_text(size = 12,color="black",face='bold',hjust=0.5,vjust=0.5),
        legend.text = element_text(size = 10,color="black",face='bold'),
        axis.text.y = element_text(size = 10,color="black",face='bold'),
        axis.title.x = element_text(size = 10,color="black",face='bold'),
        axis.title.y = element_text(size = 10,color="black",face='bold'),
        axis.text.x = element_text(size = 10,color="black",face='bold',angle = 45 , hjust=0.5,vjust=0.5) ,
        axis.line = element_line(size = 0.5)) 

out_name <- paste0( image_path , "/DSC2_DSG2.MutVsWild.IM.pdf" )
ggsave(file=out_name,plot=plot,width=width,height=5)

out_name <- paste0( image_path , "/DSC2_DSG2.MutVsWild.IM.tsv" )
write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)